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ABSTRACT 


It was the objective of this investigation to study the subsonic flow over a missile 
forebody. To this end the NASA Ames developed OVERFLOW, three dimensional, 
Navier-Stokes Code was applied and detailed results were obtained for the missile 
forebody. The flow conditions were set at a Mach Number of 0.3 and Reynolds Number 
of two million for each angle of attack at 0, 2, 6, 10 and 14 degrees. These viscous flow 
solutions were compared with three inviscid flow solutions, namely Slender Body Theory, 
von JCarman's Method and the NASA Ames Panel Code called PMARC. As expected, the 
OVERFLOW and PMARC solutions were in good agreement at zero and anaU angles of 
attack. As the angle of attack was increased, the OVERFLOW solution showed the 
development of the vortical flow separation on the leeward side of the body. 





TABLE OF CONTENTS 


I. INTRODUCTION.1 

II. SLENDER BODY METHOD FOR APPROXIMATING FLOW PAST 
BODIES OF REVOLUTION AT ZERO DEGREE ANGLE OF ATTACK... 3 

A. THEORY.3 

B. COMPUTER IMPLEMENTATION OF THE METHOD. 4 

C. MAPLE PROCEDURE VERIFICATION.5 

D. SUMMARY.7 

II. VON KARMAN'S METHOD FOR APPROXIMATING FLOW PAST 

BODIES OF REVOLUTION.9 

A. THEORY.9 

B. COMPUTER IMPLEMENTATION OF THE METHOD.11 

C. ERROR NORM ANALYSIS.12 

D. SUMMARY.14 

IV. PANEL METHOD FOR APPROXIMATING FLOW PAST 

BODIES OF REVOLUTION.15 

A. THEORY.15 

B. MODELING WAKES FROM BODIES OF REVOLUTION..17 

C. VERinCATION OF PMARC SPHERICAL BODY MODEL.18 

D. VERIFICATION OF PMARC SPINDLE BODY MODEL AT 

ZERO DEGREE ANGLE OF ATTACK.19 

E. VERIFICATION OF PMARC SPINDLE BODY MODEL AT 

nVE DEGREE ANGLE OF ATTACK.21 

F. SUMMARY.23 

V. A CFD FLOW MODEL OVER A TYPICAL MISSILE FOREBODY.25 

A. A BRIEF THEORY FOR THE OVERFLOW CODE.25 

1. Description of Basic Fluid Dynamic Equations.25 

2. Description of the Major Mathematical Operations.26 

B. OVERFLOW MISSILE FOREBODY MODEL DESCRIPTION..... 28 

vii 

























C. OVERFLOW MISSILE FOREBODY FLOW SOLUTION 

CONVERGENCE STUDY.29 

D. PMARC MISSILE FOREBODY MODEL DESCRIPTION.29 

E. COMPUTING TIMES FOR OVERFLOW AND PMARC.34 

F. OVERFLOW AND PMARC Cp DISTRIBUTION 

COMPARISONS AT VARIOUS ANGLES OF ATTACK.34 

G. SUMMARY.47 

VI. CONCLUSIONS AND RECOMMENDATIONS.49 

APPENDDC A. MAPLE PROCEDURE TO CALCULATE Cp VALUES ON A 

SPINDLE.51 

APPENDIX B. FORTRAN LISTING OF THE SPHERE CODE.61 

APPENDIX C. SPROUT FILE FROM THE SPHERE CODE.63 

APPENDIX D. ERR.OUT FILE FROM THE SPHERE CODE.71 

APPENDIX E. PMARC INPUT FILE FOR A SPHERE.73 

APPENDIX F. PMARC INPUT FILE FOR A SPINDLE.75 

APPENDK G. OVERFLOW INPUT FILE FOR A MISSILE FOREBODY.77 

APPENDIX H. PMARC INPUT FILE FOR A MISSILE FOREBODY.79 

LIST OF REFERENCES.115 


viii 


INITIAL DISTRIBUTION LIST. 


117 



















ACKNOWLEDGEMENTS 


I wish to give thanks to a number of individuals who have supported my education 
and our research goals. First, I would like to thank the Lord for giving me the strength and 
courage to undertake these enterprises in the waning years of my career. Second, to Mr. 
Milt Burford, Dr. Jim DeSanti and many others at the China Lake Naval Air Warfare 
Center, I want to give thanks for allowing me the opportunity to pursue this education. 
Third, to Professor Max Platzer and Professor Garth Hobson, I give grateful thanks. Both 
of these superb instructors have unselfishly and patiently guided my research through the 
fascinating field of computational fluid dynamics. Fourth, I could not have performed this 
endeavor without the love and cooperation of my wife Linda, our daughter Suzanne and 
her fiancd Brian, my mother Blanche and my grandmother Harriet. They have cared for 
and comforted me throughout this journey. Finally, I dedicate this thesis to my father Eric, 
who past away before seeing the completion of this venture. 


IX 


I. INTRODUCTION 


Over the years, a number of methods have been developed to provide the missile 
aerodynamicist with tools to predict the aerodynamic characteristics of various missile 
configurations. They range from purely empirical methods and simple slender body theory 
to panel methods and Navier-Stokes Methods. A recent review of these methods can be 
found in Reference 1. This review reveals the fact, that the prediction of strong viscous 
and separated flow effects, presents a great challenge. 

Modem missiles tend to have small lifting surfaces. Therefore, the contributions of 
the body to the forces and moments, and the understanding of the wing-body interference 
becomes very important. In addition, the missiles are launched at relatively high angles of 
attack, further exacerbating the problem, due to the vortex shedding from the missile 
body. As a result, there is a growing need to compute the detailed flow separation and the 
vortical flow features, which occur over typical modem missile configurations, using the 
full viscous flow equations. 

At NASA-Ames a considerable amount of work was done to study the viscous 
flow surfaces over blunt-nose bodies of revolution. This work used the NASA-Ames 
developed OVERFLOW Code for the solution of the thin-layer Navier-Stokes equations. 
References 2 through 5 document this effort. Recently Reference 6 extended this work, by 
computing of the flow field about a complete missile at high angle of attack. The results of 
Reference 6 were compared to experiments, described in Reference 7. 

As pointed out, for example in Reference 8, vortices are shed from the forward 
body and the canards. These vortices affect the pressure distribution and the forces and 
moments on the missile. Therefore, in this thesis the missile forebody of Reference 6 is 
investigated in greater detail. 

Starting with the solution for zero angle of attack, the angle of attack is gradually 
increased to 14 degrees. This allows the comparison of the Navier-Stokes solutions with 
the simpler and faster methods for the more benign flows at small angles of attack, where 
inviscid methods can be expected to give reasonable results. This comparison allows 
assessing the range of validity of the simpler methods. Slender Body Theory, von 
Karman’s Method and NASA-Ames panel code PMARC were chosen to compare with the 
Navier-Stokes calculations from the OVERFLOW Code. These methods are briefly 
described in Chapters 2 through 4. Chapter 5 presents the major results of this thesis, 
while Chapter 6 gives the conclusions and recommendations for future work. 
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n. SLENDER BODY METHOD FOR APPROXIMATING FLOW PAST BODIES 
OF REVOLUTION AT ZERO DEGREE ANGLE OF ATTACK 

A. THEORY 

The Slender Body Method is based on the assumptions of incompressible, inviscid 
and irrotational flow. The governing equation is Laplace’s equation for the potential 
(t)(x,r,0) [Reference 9], 


^ ^ 1^ 1 

dx^ dr^ r dr 30 ^ 


( 2 . 1 ) 


Equation 2.1 is in cylindrical coordinates (x,r,0). The body's angle of attack is 
assumed to be zero, and the body's longitudinal axis is assumed to coincide the x axis. 

Source superposition along the x axis, then produces the well-known equation for 
the disturbance potential. 


45iJr 


0 r 


f(4)4 

4 ( 4 -.)^]” 


( 2 . 2 ) 


(U) is the free stream speed. (^) is the source location on the body’s longitudinal axis. The 
body extends from 0 to /. (F) is the change of the body's cross sectional area along the 
longitudinal axis. It is given by. 


F{x) = 2nR{x)^^ 
dx 


The body radius is given by. 


r = F(x) 


(2.3) 


(2.4) 
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Depending on the complexity of the surface function, the integrals in Equations 2.5 
and 2.6 can be difficult to integrate analytically. Therefore, the MAPLE Program was used 
to symbolically integrate Equations 2.5 and 2.6 and solve Equation 2.7. 

B. COMPUTER IMPLEMENTATION OF THE METHOD 

MAPLE is a general purpose, commercially available, computer code that can 
perform symbolic mathematics [Reference 11]. The code may also be used for numerical 
mathematics. To demonstrate this program's usefulness in potential flow theory. Appendix 
D lists the implementation of Equations 2.5, 2.6 and 2.7 in Sections 1 through 3 
respectively. The MAPLE Procedure, given in these Sections, models the flow over a 
spindle body of revolution as shown in Figure 1. The spindle is at zero degree angle of 
attack. 
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Figure 1. Spindle Body of Revolution Surface Function and Shape. 


The geometry is non-dimensional, so that the length is one and maximum diameter 

(tau) is 16% of the body's length. The maximum diameter occurs at the middle of the 
body. Velocity components and values for 10 points on the spindle's surface are given 

in Section 3 of Appendix A. 

C. MAPLE PROCEDURE VERIFICATION 

The parabolic spindle shown in Figure 1 is described by, 

R{x) = 2tau{x-x^) (2.8) 

Its first derivative becomes, 

^^^ = 2tau{l-2x) ( 2 . 9 ) 

dx 

The absolute magnitude of Equation 2.9 will always be less than one, since the maximum 
diameter is 0.16. This satisfies the slender body approximation. 

To verify the spindle results in Section 3 of Appendix A, Reference 10 gives an 
approximation of Equation 2.5 of the form. 
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In this equation (F ) is Equation 2.3, while (F ), (F") and (F'”) represent second, third 
and fourth derivatives of the body's cross sectional area respectively. Section 4 of 
Appendix A determines velocity values for Equation 2.10. These are compared to values 
derived from Equation 2.5, using the surface function for the spindle in Figure 1. Figure 2 
shows the comparison between Equations 2.5 and 2.10 for 48 surface points. The 
agreement is seen to be excellent. Both solutions exhibit the well-known singular behavior 
of slender body theory at the nose and tail of the body. 



Figure 2. Comparison of Equation 2.5 and 2.10. 
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D. SUMMARY 


This singularity behavior is the result of the integration of Equations 2.5 and 2.6 in 
Sections 1 and 2 of Appendix A respectively. In both cases integration yields terms 
containing the logarithmic expression, 

ln(2Vx^+r^ -2x) (2.11) 

At the spindle's leading or trailing tip the radius approaches zero. At the tips the x 
coordinate approaches either zero or one respectively. Under these conditions the 
logarithmic expression becomes indeterminate. As depicted in Figure 2, the velocity 

components do not have a finite value at the tips of the body. Expression 2.11 also 
appears in Equation 2.7 (see Section 3 of Appendix A), indicating that the at the tips 

will not be the expected finite value of one. Another indication of this indeterminate 

behavior is given in Sections 5 and 6 of Appendix A. 

In Section 5 as r and x approach zero, the value peaks at approximately 0.89 

The values decrease and do not approach the expected value of one. In Section 6 as r 
approach's zero and x approach's one, the value oscillates about 0.64, then suddenly 

decreases at the x coordinate of 0.99927. At this point the MAPLE Program reports a 
singularity in one of the logarithmic expressions occurring in the solution of Equation 2.7. 

This implies that the Slender Body Method is inherently inaccurate near the leading 
and trailing tips of a body of revolution. The method can be used to confirm 
computational fluid dynamic (CFD) models in the middle portion of the body, provided 
that the surface function is relatively simple and continuous. When surface functions can 
not be easily formulated for complex shapes, the Panel Method might prove more useful. 
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in. VON KARMAN'S METHOD FOR APPROXIMATING FLOW PAST 

BODIES OF REVOLUTION 


A. THEORY 


von Karman's Method assumes a potential flow about a body of revolution 
[Reference 12]. The flow is considered incompressible, inviscid and irrotational. The 
governing equation is Laplace's equation, 


ay lay i 

dz^ dr^ r dr r' 80' 


(3-1) 


This represents the flow within a body of revolution in cylindrical coordinates (z,r,6) 
[Reference 13]. The solution to this equation may be formed by superposition of 
elementary flows. These are represented by stream functions (\|r). The stream function for 
uniform flow of speed (U) z direction is given by, 

= (3.2) 


To obtain a solution for flow over a body of revolution at zero angle of attack (non-lifdng 
flow), whose axis coincides with the z axis, it is sufficient to distribute sources along the 
body axis. The stream function of a source located in the origin of the coordinate system 
is. 


V = - 


mz 


(3.3) 


The change in stream function ) at any point in the flow from sources along a 
small interval (dq) there is. 



qj{z-q)dq 

^lr^+{z-qf 


(3.4) 
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The interval is measured along the body's longitudinal axis, is the source strength per 

unit length. Integration of Equation 3.4 along this axis yields the induced stream function 
at any point in the flow. This function is created by sources on the longitudinal axis. The 
effect of all the source strengths and uniform flow are summed to give the total stream 
function at any point in the flow. From Reference 13 the summation can be represented as, 


>|((r,z) = it/r=+X 


+(z-Zj^,y -yjr^ +iz-z'jy qj (3.5) 


The (zj+i) and (z^) coordinates define the location of the source of constant strength (qj) 

on the longitudinal axis. The body of revolution is defined by Equation 2.4. Realizing that 
the body surface represents the stagnation streamline \j/ = 0, equating the right hand side 
of Equation 3.5 to zero provides a means to determine the unknown source strength qj. 

For finite length bodies, the total sum of the source strengths must be zero. This 
additional condition must be appended to the system given by Equation 3.5. The system 
therefore becomes. 


{const}=[dd]{q} 


(3.6) 


The vector {const} is only a function of the known body radius. The coefficient matrix 
[dd] is a function of the known body geometry and source locations, while the vector (q) 
is the unknown source strengths. The resulting system of algebraic equations may now be 
solved by common equation solves , such as Cramer's Rule, Gaussian Elimination, etc. 

Once the source strengths are known on the body's surface, the flow velocity 
components (u^, u^) can be defined through the first derivative of Equation 3.5. From 
Reference 13 the first derivative with respect to z yields the radial component u^,. The 
longitudinal component is obtained from the first derivative with respect to r. Using 
Bernoulli's Equation, the non-dimensional coefficient of pressure (C^) is obtained as. 



(3.7) 


The Cp is the ultimate result of this Method. Since the is a function of the bod/s 
geometry, the value will differ at various points in the outer flow field and on the body's 
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surface. This condition makes the parameter an excellent metric to make comparisons 
of results from various analytic and computer code techniques. 

B. COMPUTER IMPLEMENTATION OF THE METHOD 


Reference 13 lists a computer code that implements von Karman's Method. The 

code was called SPHERE. A FORTRAN listing of the code is given in Appendix B. This 
code predicts values about a sphere of radius (a). Since the flow field is symmetrical, 

the code only calculates values at surface points on the upper half of the sphere. The 

free stream velocity and the radius of the sphere were assumed to be one. SPHERE used 
Cramer’s Rule to solve the system of equations. 

Appendix C and D contain listings of the two output files created by SPHERE. 

The file called SPH.OUT in Appendix C gives the surface point location (r^, z,.), 
calculated (C^) and exact {CpCx) for various numbers of surface points on the sphere. 
The exact values are known to be given by, 


CpCx = 1 


IfnX 

4UJ 


(3.8) 


The derivation of Equation 3.8 may he foimd in Reference 13. Calculated and exact 

values were determined for nine surface points as a check of the code. Figure 3 shows 
that the calculated values plot identically on the exact Cp values for the nine points. 

However, when the number of surface points is increased, the calculated Cp values 

diverge from the exact. This is expected because the matrix becomes increasingly ill- 
conditioned. 
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Figure 3. Calculated and Exact Cp Distributions for 9 Surface Points on a Sphere. 
C. ERROR NORM ANALYSIS 

Figures 4 and 5 show how the calculated Cp distributions diverge from the exact 

distribution, as the number of surface points increases. This divergent behavior comes 
from the ill-conditioning of the coefficient matrix [dd] in Equation 2.6. The growth of this 
ill-conditioning can be defined as, 

{error}=[dd]{q} - {const} (3.9) 

error norm = ^'^{errory (3.10) 

Appendix D lists the ERR.OUT file from the SPHERE Code. This listing gives the 
results of Equation 3.10 for various numbers of surface points. Examination of this data 
and the data's plot in Figure 6 clearly shows the error growth, as the number of surface 
points increases When the number of surface points exceeds 39, the matrix [dd] becomes 
so ill-conditioned that the code cannot obtain the vector {q) in Equation 3.6. 
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Figure 4. Calculated and Exact Cp Distributions for 33 Surface Points on a Sphere. 
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Figure 6. Error Norm verses Number of Surface Points for a Spherical Body of 

Revolution. 


D. SUMMARY 

As shown by von Karman, this method of representing incompressible, inviscid 
flow over a sphere can be extended to the flow over arbitrary bodies of revolution at zero 
and small angles of attack. However, the method is limited to slender bodies, because 
source or doublet distributions are used only along the axis. In contrast to the slender 
body theory described in Chapter H., it has the advantage of describing the flow over 
blunt-nosed bodies. 

A more general method that eliminates the restriction to slender bodies is afforded 
by surface distributions of sources and doublets. Such a method is described in the next 
chapter. 
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IV. PANEL METHOD FOR APPROXIMATING FLOW PAST BODIES OF 

REVOLUTION 


A. THEORY 


The Panel Method is also based on potential flow theory and generally parallels the 
development of the Slender Body and von Karman Methods. Reference 14 gives a detailed 
mathematical description of the Panel Method. This reference also describes how the 
Method is implemented in the computer code called PMARC (Panel Method Ames 
Research Center). PMARC assumes an incompressible and irrotational flow. The code 
does allow a viscous boundary layer on the body, but this option was not used. All of the 
PMARC results in this thesis were based on inviscid flows. The method solves the three 
dimensional form of Laplace's equation, 

V2(t) = 0 (4.1) 


As shown in Reference 14, application of Green's Theorem leads to the following 
integral equation. 



(4.2) 


(|i) and (|J,„,) are unknown doublet strengths, on the body's surface element (dS) and on 
the wake (TV) respectively. The doublet strength (|Xp) is zero, if the arbitrary point (P) of 
interest is off the surface. When the point is on the surface, the strength is 27c. (n) is the 
outward surface normal on the body's surface, while (r) is the vector between P and dS. 
n and r are known from the body's geometry. The source strength (a) is given by, 

(4.3) 


On the body's surface the normal velocity is usually zero or may be prescribed. The 
free stream velocity vector (1^) is known. Therefore, the source strength is known in 
Equation 4.2. 
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Equation 4.2 may be solved numerically by discretizing the body's surface into 
panels. These panels represent areas of source and doublet strengths, that when known, 
will define the flow field over the body's surface. This marks the only difference between 
the Panel Method and those given in Chapters n and HI. Previous methods located source 
strengths inside the body on its' longitudinal axis. In PMARC the source and doublet 
strengths are assumed to be constant over a panel. 

The discretization process results in a system of algebraic equations of the form, 

Ns \ / \ 

(4-4) 

k=l it=l /=! 

The summation terms in Equation 4.4 are analogous to the integrals of Equation 4.2. (B) 
and (O become respectively. 


k ' 

(4.5) 


(4.6) 


The first and second terms in Equation 4.4 represent the placement of surface panels on 
the body. The third term accounts for the wake panels trailing downstream from the body. 
In PMARC this equation is solved by an iterative technique for very large matrices. The 
solution occurs in a time stepping fashion, so that either steady or unsteady flows may be 
calculated. The wake doublet strengths become a function of the surface doublet strengths 
or, 


(4.7) 

This means that the PMARC time stepping algoritiim effectively substitutes Equation 4.7 
into Equation 4.4 to yield the surface doublet strengths directly. 

These doublet strengths are differenced to obtain the tangential velocity 
components on the surface of each panel. The velocity components are then used in. 
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(4.8) 



^87t'l 

rp(r)-4(i-l)^ 

v2 


V At j 


Equation 4.8 defines the pressure coefficient (C^^) at each panel (k). (V^) is the flow 

velocity magnitude at each panel. The third term in Equation 4.8 accounts for pressure 
changes when the flow is unsteady. This term is zero for steady flows, so that Equation 
4.8 becomes identical to Equations 2.7 and 3.7. Because the theory and PMARC allow the 
effects of wakes, a problem arises as to how to attach wake panels to bodies of revolution 
to simulate the effect of flow separation. 

B. MODELING WAKES FROM BODIES OF REVOLUTION 

PMARC allows wakes to be excluded from the model. When this occurs PMARC 
will calculate Cp values, but can not calculate circulation about the body. In this case total 

aerodynamic forces, such as lift, drag, etc. are not determined. Alternatively, when wakes 
are included in the model, they must be defined as tither flexible or rigid. Flexible wakes 
require that the separation line be specified in the input file. Rigid wakes require initial 
wake panels be defined in regions where the flow is expected to separate. 

Reference 14 strongly implies that PMARC's wakes should only be modeled when 
attached to traditionally defined lifting surfaces, that have sharp, linear trailing edges. This 
allows a precise definition of a separation line or initial wake panel layout along the edge. 
The code may then properly apply the Kutta condition at the edge. When wakes are 
attached to bodies of revolution, without lifting surfaces, the results appear to be 
questionable, depending on whether the body has a tapered, hemispherical or flat base 
[Reference 15]. A literature search of References 16 through 19 provided no guidance into 
acceptable wake modeling for bodies of revolution. Reference 18 claimed to have modeled 
wakes on a body with no lifting surfaces, but provided no details of the PMARC input file. 

In view of this situation the author did not include wakes in any of the PMARC 
flow models in the thesis. All flow models were considered to be steady. This meant that 
the third terms of both Equations 4.4 and 4.8 had a value of zero. To evaluate PMARC’s 
ability to compute fully attached flows the code was first applied to the sphere and the 
spindle, so that the results could be compared with the Slender Body and von Kaiman's 
Methods. 
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c. 


VERIFICATION OF PMARC SPHERICAL BODY MODEL 


A spherical body was constructed for the PMARC Code. Appendix E lists the 
input file for the sphere. In the file the name lists &BINP8, &BINP8A and &BINP8B 
were set to zero in the file. These settings created steady flow. The parameters in the 
«&:WAKE1 command were also set to zero. This eliminated the wakes from the model. 
Figure 7 shows the surface paneling layout. 



Figure 7. Surface Panel Layout for a Spherical Body. 


The panels were positioned on the sphere to yield nine values that 
corresponded to identical radial and longitudinal locations of the exact Cp values given in 

Reference 13. AU PMARC results shown in this thesis were generated by version 12.21, 4 

March 1994, of the code. Figure 8 shows excellent comparison between PMARC and 
exact Cp values. 
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Figure 8. PMARC and Exact Distributions for nine Surface Points on a Sphere 

without wakes. 

D. VERIFICATION OF PMARC SPINDLE BODY MODEL AT 
ZERO DEGREE ANGLE OF ATTACK 

This trend was also confirmed when the spindle body was modeled in steady flow, 
without wake panels and at zero degree angle of attack. A typical PMARC input file for 
the spindle is given in Appendix F. Figure 9 displays the PMARC spindle model, using the 
surface function given in Figure 1. In Figure 9 a positive rotation about the y-axis 

produces a nose up positive angle of attack, using the right hand rule. The geometry was 
non-dimensionalized similar to Figure 1. values were calculated and compared to 

similar values fi'om the Slender Body Method, described in Chapter n. Figure 10 shows 
that PMARC results compare favorably. 
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Figure 9. Surface Panel Layout for a Spindle Body. 



Figure 10. PMARC and Slender Body Distributions for 48 Surface Points on a 
Spindle at Zero Degree Angle of Attack mtiiout Wakes. 
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E. 


VERIFICATION OF PMARC SPINDLE BODY MODEL AT 
FIVE DEGREE ANGLE OF ATTACK 


The PMARC spindle model was exercised at both zero and five degree angles of 
attack. This was done to verify the model's physical behavior from a point of view. 

Figure 11, 12 and 13 compare distributions on longitudinal lines along the body's 

suction or top surface, middle or mid plane surface and pressure or bottom surface 
respectively. These longitudinal lines can, alternatively, be defined in Figure 9 as lying in 
the positive x, z plane for the top surface, the positive x, y plane for the mid plane surface 

and the positive y, negative z plane for the bottom surface. As expected for a symmetrical 
spindle, at zero degree angle of attack the distributions are the same longitudinally, 

regardless of surface orientation. 

When the spindle is elevated to an attack angle of five degrees, suction or negative 
pressures develop on the forward top surface in Figure 11. As expected, the flow 
accelerates over the forward top surface, creating negative pressures. On the aft top 
surface the flow experiences an adverse pressure gradient, slowing the flow and creating 
positive pressure. 



Figure 11. PMARC Distributions for 48 Surface Points on a Spindle at Zero and 
Five Degrees' Angles of Attack on the Top Surface without Wakes. 
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Figure 12. PMARC Distributions for 48 Surface Points on a Spindle at 2^ro and 
Five Degrees' Angles of Attack on the Mid Plane Surface without Wakes. 



Figure 13. PMARC Distributions for 48 Surface Pmnts on a Spindle at Zero and 
Five Degrees' Angles of Attack on the Bottom Surface without Wakes. 
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Figure 12 indicates a relatively constant difference between values at attack 

angles of zero and five degrees along the spindle's mid plane surface. At zero degree a 
fluid particle experiences only body fixed, x and y velocity components, as it travels along 
the body's mid plane surface [see Figure 9]. The particle travels circumferentially and 
longitudinally along the body's side, when the body is rotated to five degrees. In this case 

the particle's velocity vector contains all three velocity components in the x, y and z 
directions. The total velocity is different for either case, producing different 

distributions. 

As expected, in Figure 13 the flow slows on the forward bottom surface, creating 
positive pressures. On the aft bottom surface the flow experiences a favorable pressure 
gradient, creating negative pressures. 

F. SUMMARY 

The Panel Method, as implemented in the PMARC Code, predicts distributions 

for bodies of revolution that compare very well with classical potential theories. The code 
predicts the expected physical behavior of bodies at small angles of attack. Only steady, 
incompressible, inviscid and irrotational conditions should be considered in the flow. More 
investigation is needed to understand PMARC's predictions very near the tips of conical 
shaped bodies. Guide lines should be developed for modeling wakes on bodies of 
revolution. 

Even with these slight deficiencies, PMARC is a suitable tool to confirm pressure 
fields about CFD models, as long as viscous effects are small. Before explaining how 
PMARC may be used to verify these models, a brief discussion of the CFD computer code 
used in this study will be given. 
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V. A CFD FLOW MODEL OVER A TYPICAL MISSILE FOREBODY 
A. A BRIEF THEORY FOR THE OVERFLOW CODE 
1. Description of Basic Fluid Dynamic Equations 

Due to the experience discussed in Chapter I, the Navier-Stokes flow solver called 
OVERFLOW was chosen for this study. OVERFLOW is currently in development at the 
NAS A-Ames Research Center. 

From Reference 20 the governing equations may be written in the three Cartesian 
X, y and z coordinates. In the vector form the equation for conservation of the flow mass 
becomes, 

^+V»py = 0 (5.1) 

at 

The conservation of flow momentum in compact form is, 

pf = V.C, (5.2) 

Three momentum equations are obtained, when Equation 5.2 is expanded in each 

Cartesian direction. The conservation of energy equation in the flow may be written as, 

p^=v.(tvr)+o,|i (5.3) 

Dt dXj- 

In these five equations (p), (t), (V), (e^), (k) and (T) are respectively density, time, the 

velocity vector, internal energy, thermal conductivity and temperature. The fluid stresses 
(c,y) are split into pressure and viscous terms as. 
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(5.4) 


<y,7= A+4 


I 

^dxj dXi j 




The first term in Equation 5.4 is the normal pressure (p) and (5,y) is the Kronecker Delta 


Function. This term represents the normal stress components on the fluid element. In the 


second term (p.) is the molecular viscosity and 


dX; 




is the strain in the fluid. In the 


third term (^) is the bulk viscosity and is usually assumed equal to [ p | [Reference 


21]. When this assumption is made, the last two terms in Equation 5.4 may be combined 
to form the viscous or shearing stress components on the fluid element. 

In OVERFLOW the molecular viscosity in the flow field is approximated by 
Sutherland's Formula as, 


p rrf'Yr,+5>[ 
vio [tJ Vt+s) 


(5.5) 


(p^) and (7^) are the reference flow viscosity and temperature respectively. The parameter 
(5) is assumed to be equal to 194 degrees Rankine for standard air. The turbulent 
viscosity is determined by the Baldwin and Lomax Turbulence Model in the boundary 
layer. This turbulent viscosity is added to Equation 5.5. Other turbulence models are 
available in the code [Reference 20]. 

Equations' 5.1, 5.2 and 5.3 allow modeling of compressible, viscous and vortical 
flows, removing the basic assumptions underlying the previous potential flow theory 
methods. Various mathematical operations are performed on Equations 5.1, 5.2 and 5.3. 
This effectively places them in a finite difference form suitable for a computer solution. 


2. Description of the Major Mathematical Operations 

Equations' 5.1, 5.2 and 5.3 are placed in conservative vector form, in order to 
maintain the accuracy of the solution. The equations become. 
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(5.6) 


dQ dF dG dH ^ 
dt dx dy dz 

Q = [p,pu,pv,pw,ef 

F = [pu,pu^ +p-x„,puv-x^,puyv-x^Xe+p)u-ux^-vx^-wx^ 

G = [pv, pvu-x^,pv^+p-Xyy,pvw- Xy, ,{e + p)v-ux^ - vx^ - wx^ + f 

H = [pw,pwu-x^,pwv-x^,pw^ + p-x^,{e+p)w-ux^ - vx,, - wx^ + qj 

(u ), (v) and (w) are the velocity components. (x,y, i,j = x,y,z) are the viscous stress 

components. These stresses can be formed from the last two terms in Equation 5.4. (e) is 
the total energy and (q^, i = x,y,z) are the heat transfer terms. Equation 5.6 is non- 

dimensionalized, so that the different flows may easily be compared, using parameters, 
such as Mach, Reynolds and Prandtl Numbers. Solutions from non-dimensionalized 
equations may also be applied to a flow over different size bodies of the same shape. The 
equations are transformed from the physical domain into computational domain. This 
transformation is done, so that the equations may eventually be differenced on a uniformly 
spaced grid [Reference 22]. Equation 5.6 is further simplified by the thin layer 
approximation, which neglects the viscous terms over the body in the streamwise 
direction.. The equations can now be cast into finite difference form. This leads to a 
system of algebraic equations, which are solved implicitly in OVERFLOW. 

The system in a very simple form may be written as, 

[/h5]{AQ} = {r;M} (5.7) 

This is called the delta form [Reference 4]. The matrix [/hs] and vector {rfe} are 
functions of known flow quantities, such as density, pressure, etc. The vector {A0} is the 
unknown change in the dependent variables between the current and next time step. For 
steady flows, when the change in these variables approaches zero, the solution is 
considered to have converged to the best approximation of the original Equations 5.1 
through 5.3. 
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B. OVERFLOW MISSILE FOREBODY MODEL DESCRIPTION 


A flow field about a missile forebody was constructed and exercised in 
OVERFLOW to obtain an acceptable converged solution. Figure 14 shows the missfle 
forebody's surface, but for clarity it does not show the flow grid extending above the 
surface. The missile is oriented in an Z type configuration to the flow. The axis system is 
body fixed with the same angle of attack definition given in Figure 9. All of the dimensions 
of the grid have been non-dimensionalized with respect to the missile's maximum diameter. 
The maximum diameter is one, while the missile's length is approximately 17.5 diameters. 
The flow grid extends 39 diameters forward of the nose. The grid extends over the missile 
body's surface at the aft end to 75 diameters and for computational efficiency only one half 
of the flow field is modeled. The flow is assumed symmetric about the x, z plane. There is 
no grid in the negative y direction. There is no grid aft of the missile body's rear base. 

The grid over the missile forebody's surface is composed of six sub-grids. These 
define the in-flow over the missile forebody's nose, the body's mid section, the top and 
bottom canards, and the out-flow over the body's rear section. The sixth grid is an outer 
flow grid and covers the mid section. These six grids provide 874,125 total grid points to 
describe the flow field over the missile's forebody. 

All these grids overlap each other to varying extents. In the overlapped regions the 
chimera scheme is used to provide interpolation and grid blanking information [Reference 
23]. This information is used to allow OVERFLOW to interpolate the flow variables from 
one grid to the next. The six grids and the interpolation information are provided to the 
code in two large files called GRID.IN and FORT.2 respectively. The files were from 
Reference 24. In addition to these files OVERFLOW requires an input file that describes 
such things as the number of time steps desired, angle of attack, Mach and Reynolds 
numbers, time step size, boundary conditions, turbulence model, etc. Appendix G lists a 
typical OVERFLOW input file for the missile forebody. A description of the various input 
parameters can be found in Reference 20. All OVERFLOW results shown in this thesis 
were derived from version L6ao, 8 February 1994, of the code. The forebody model was 
initially exercised at a benign flow condition of zero angle of attack, in order to check the 
flow physics and solution convergence. 
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c. 


OVERFLOW MISSILE FOREBODY FLOW SOLUTION 
CONVERGENCE STUDY 


The initial flow field was set for a free stream Mach Number of 0.3, a Reynolds 
Number of 2.0 milli on and zero degree angle of attack. Figures 15 through 20 show the 
plots of the norm of the vector {Aj2} in Equation 5.7 at each iteration. The norm is 
the square root of the sum of squares of the vector's elements. As the norm approaches 
zero, the change in the dependent variable also approaches zero, creating a converged 
solution. The plots in these figures indicate an acceptable converged solution, after 4000 
time steps. 

D. PMARC MISSILE FOREBODY MODEL DESCRIPTION 

A PMARC flow field of the missile forebody was constructed, using 4836 surface 
panels. As discussed in Chapter IV, wake panels and boundary layers were not included. 
The geometry of the PMARC missile forebody in Figure 21 was very similar to the 
OVERFLOW forebody in Figure 14. Some slight differences existed between the canards 
of the two forebodes. The PMARC model uses the same body fixed coordinate system 
and angle of attack definition in Figures 9 and 14. Appendix H lists a typical input file for 
the PMARC missile forebody. A description of the various input parameters can be found 

in Reference 14. The PMARC model was exercised at various angles of attack and 
compared to distributions derived from the OVERFLOW forebody model. 

Figure 22 compares typical distributions from both the OVERFLOW and 
PMARC models. This figure indicates that most of the changes occur in the nose 

region of the missile. This region extends approximately three to five calibers aft the nose, 

along the x axis. Figure 23 shows the pressure distribution over the first three calibers and 
the missile body profile. Similar Cp distribution comparisons were made in this region at 

various angles of attack, in order to demonstrate where the two models begin to differ. 
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Figure 14. OVERFLOW Missile Forebody Surface. 
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Figure 17. Norm for the Outer Flow Grid. 



Figure 18. Norm for the Rear Section Grid. 
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Figure 19. Norm for the Top Canard Grid. 
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Figure 20. Norm for the Bottom Canard Grid. 
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E. COMPUTING TIMES FOR OVERFLOW AND PMARC 


Typically OVERFLOW calculations required approximately 10 to 12 cpu days on 
the Naval Postgraduate School, Cray Y-MP/EL computer for each converged solution. A 
typical PMARC solution required 1 to 2 cpu days on a Silicon Graphics workstation. The 
Cray mass storage disk space was accessed across the network for the PMARC runs, 
since local workstation space was very limited. A PMARC solution could have been 
obtained much faster, if more local disk space were available. 

F. OVERFLOW AND PMARC Cp DISTRIBUTION 

COMPARISONS AT VARIOUS ANGLES OF ATTACK 

Figures 25 through 38 present distributions for the top , mid plane and bottom 

surfaces of the nose regions of the two missile forebody models. The models were 
exercised at 0,2, 6,10 and 14 degree angles of attack. The scale used in these figures was 
the same as used in Figure 23. Figures 39 and 53 are enlargements of figures 25 through 
38 in the nose region. These figures show how OVERFLOW and PMARC results diverge 
as the angle of attack increases. 

At zero and two degree angles of attack in Figures 24 through 29 and Figures 39 
through 44, there are very good agreement between the two models. At six degree angle 
of attack in Figure 30 through 32 and Figures 45 through 47, some differences begin to 
appear along the mid plane surface. The differences in the aft portion of the bottom 
surface are significant in Figure 32. 

Figures 33 through 35 and 48 through 50 show distributions, when the models 

are at 10 degree angle of attack. Both models are in excellent agreement within the first 
caliber along the x axis. However, Figures 33 through 35 show differences between both 
models in the aft portion of the nose region. 

Figures 36 through 38 and 51 through 53 present the results of the two models at 
14 degree angle of attack. These figures show agreement within the first caliber along the 
X axis, similar to the 10 degree case. Again there are differences between the 
OVERFLOW and PMARC models at 14 degree angle of attack in the aft portion of the 
missile's nose. 
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Figure 21. PMARC Missile Forebody Surface. 





Figure 22. OVERFLOW and PMARC Distributions for the Full Missile 

Forebody. 



Figure 23. OVERFLOW and PMARC Distributions for the Missile Nose. 
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Figure 24. OVERFLOW and PMARC Distributions for the Missile Forebody at 
Zero Degree Angle of Attack on the Top Surface. 
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Figure 25. OVERFLOW and PMARC Distributions for the Missile Forebody at 
Zero Degree Angle of Attack on the Mid Plane Surface. 



Figure 26. OVERFLOW and PMARC Distributions for the Missile Forebody at 
Zero Degree Angle of Attack on the Bottom Surface. 
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Figure 27. OVERFLOW and PMARC Distributions for the Missile Forebody at 
Two Degree Angle of Attack on the Top Surface. 


Figure 28. OVERFLOW and PMARC Distributions for the Missile Forebody at 
Two Degree Angle of Attack on the Mid Plane Surface. 


Figure 29. OVERFLOW and PMARC Distributions for the Missile Forebody at 

Two Degree Angle of Attack on the Bottom Surface. 
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Figure 30. OVERFLOW and PMARC Distributions for the Missile Forebody at 
Six Degree Angle of Attack on the Top Surface. 



Figure 31. OVERFLOW and PMARC Distributions for tiie Missile Forebody at 

Six Degree Angle of Attack on the Mid Plane Surface. 
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Figure 32. OVERFLOW and PMARC Distributions for flie Missile Forebody at 
Six Degree Angle of Attack on the Bottom Surface. 
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Figure 33. OVERFLOW and PMARC Distributions for the Missile Forebody at 
Ten Degree Angle of Attack on the Top Surface. 



Ten Degree Angle of Attack on the Mid Plane Surface. 



Ten Degree Angle of Attack on the Bottom Surface. 
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Figure 36. OVERFLOW and PMARC Distributions for the Missile Forebody at 
Fourteen Degree Angle of Attack on die Top Surface. 



Figure 37. OVERFLOW and PMARC Distributions for the Missile Forebody at 
Fourteen Degree Angle of Attack on the Mid Plane Surface. 



Figure 38. OVERFLOW and PMARC Distributions for the Missile Forebody at 

Fourteen Degree Angle of Attack on the Bottom Surface. 
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Figure 40. OVERFLOW and PMARC Distributions for the Missile Forebody at 
Zero Degree Angle of Attack on the Mid Plane Surface. 



Figure 41. OVERFLOW and PMARC Distributions for the Missile Forebody at 
Zero Degree Angle of Attack on the Bottom Surface. 
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Figure 42. OVERFLOW and PMARC Distributions for the Missile Forebody at 
Two Degree Angle of Attack on the Top Surface. 



Figure 43. OVERFLOW and PMARC Distributions for the Missile Forebody at 
Two Degree Angle of Attack on the Mid Plane Surface. 



Figure 44. OVERFLOW and PMARC Distributions for the Missile Forebody at 
Two Degree Angle of Attack on tiie Bottom Surface. 
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Figure 45. OVERFLOW and PMARC Distributions for the Missile Forebody at 
Six Degree Angle of Attack on the Top Surface. 
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Figure 46. OVERFLOW and PMARC Distributions for the Missile Forebody at 
Six Degree Angle of Attack on the Mid Plane Surface. 



Figure 47. OVERFLOW and PMARC Distributions for the Missile Forebody at 

Six Degree Angle of Attack on the Bottom Surface. 
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Figure 48. OVERFLOW and PMARC Distributions for the Missile Forebody at 
Ten Degree Angle of Attack on the Top Surface. 



Figure 49. OVERFLOW and PMARC Distributions for the Missile Forebody at 
Ten Degree Angle of Attack on the Mid Plane Surface. 



Figure 50. OVERFLOW and PMARC Distributions for the Missile Forebody at 
Ten Degree Angle of Attack on tiie Bottom Surface. 
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Figure 51. OVERFLOW and PMARC Distributions for the Missile Forebody at 
Fourteen Degree Angle of Attack on the Top Surface. 
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Figure 52. OVERFLOW and PMARC Distributions for the Missile Forebody at 
Fourteen Degree Angle of Attack on the Mid Plane Surface. 



Figure 53. OVERELrOW and PMARC Distributions for the Missile Forebody at 
Fourteen Degree Angle of Attack on tiie Bottom Surface. 
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G. SUMMARY 


Missile forebody models were constructed for use in both the OVERFLOW and 
PMARC Codes. Based on distributions, comparisons were made between these two 

models for a low Mach number at several angles of attack. Agreement was very good at 
very low angles of attack. At six degree and higher angles of attack, differences between 
these two models began to appear. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


In this thesis the NASA-Ames developed OVERPLOW Navier-Stokes Code was 
used to obtain solutions for Mach = 0.3 and Reynolds = 2.0e+06 flow over a typical 
canard configured missile forebody. In the absence of detailed experimental flow and 
pressure data, the Navier-Stokes computed pressure distributions were compared with 
potential flow solutions, i.e.. Slender Body Theory, von Karman's Method and the NASA- 
Ames panel code PMARC. As expected, the first two methods were only of limited 
applicability, but PMARC provides good agreement for zero and small angles of attack. 
The PMARC and OVERFLOW computed pressure distributions started to deviate in fi-ont 
and close to the canards. This indicates the presence of flow separation and vortical flow 
development. Further, work is necessary to display the flow features in more detail, so that 
the underlying flow physics can be better understood. 


49 






50 



APPENDIX A. MAPLE PROCEDURE TO CALCULATE Cp VALUES ON A 

SPINDLE 


> 


> # This MAPLE Procedure calculates the Cp values about a spindle body of 

> # revolution at zero degree angle of attack, using potential 

> # theory. The Cp values are determined in Sections 1 through 3. 

> 

> # Section 4 determines the non-dimensional longitudinal velocity of 

> # Equation 2.10 and compares it to Equation 2.5 to confirm the 

> # procedure. 

> 

> # Sections 5 and 6 investigate the indeterminate behavior of the Cp or 

> # Equation 2.7, as the spindle’s radius approaches zero and the x 

> # coordinate either approaches zero or one respectively. 

> 


> 

># SECTION 1. Calculate non-dimensional longitudinal velocity component, 
># using Equations 2.8,2.3 and 2.5. 

> 

>tau;=; 0.16; 


T 


.16 


>R:= 

> derivative Rrs 


2.0*tau*((xi)-(xi''2)); 

/?:=.320|-.320|^ 

diff(R,xi); 


derivative Jt. := .320 - .640 ^ 

>eq_2_3:= 2.0*pI*R*derlvative_R; 

eqJJ := 2.0 Jt (.3201 - 320 (.320 - .640|) 

> lntegrand_eq_2_5:= (eq_2_3 * (x-xi)) / (((x-xi)*2 + r*2 )*(3/2)); 

> Eq_2_5:= ( (1) / (4.0*pi)) * (lnt( Integrand_eq_2_5, xis0.0..1.0 )); 

EqJJ := .2500000000 ( -.2048000000 n 

(-21.0 - 6. jc %3 + 15.00a: - 3,000 + 6.3?%3j^+ 9.xr^ + 9.j?- 3.0 + %3 ~ 3. t^ %3 J^) 

/y^+.2048000000 Jt 

Jx? + T^)ln 


%1 :=ln(z Jx^ + r^ -2.x) 


%2:=a?.2.0x+1.00 + /^ 

%3 := In(2. J%2 + 2.0 - 2. x) 

> _ 

> 

> # SECTION 2. Calculate non-dimensional radial velocity component, 

> # using Equation 2.6. 
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> integrand_eq_2_6:= (eq_2_3) / {((x-xl)''2 + r''2 )^(3/2)); 

. , n(.320^-.320^")(.320-.640^) 

integrand^q 2_o := 2.0- - - — 

> Eq_2_6:= ((r) / (4.0*pl)) * {int{ integrand_eq_2_6, xi=0.0..1.0)); 

EqJjS := .2500000000 r (-.2048000000n(4.00 + 4.Qx^+ 7.0xr^■ 5.0:^+ 2.x* 

.6.xln{2.Jm+2.0-2.x)r^Jm-4.i^-1.0x + 3.i^ln{zJ^ + 2.0-2.x]J^yir^Jm] + 
.2048000000n(2.X* + /^-2. J(?+ 3.ln(2. Jx^ + - 2.x) Jx^ + 

-6.xln(2. + -2.x) -4.r*-3.xr^-3.j^)l{r^J^ + r^))(n 

%1 ;=a?-2.0x+ 1.00 + /^ 


> # SECTION 3. Calculate the Cp values, 

> # using Equation 2.7. 

> Cp:= ( -2.0*(Eq_2_5)) - (((Eq_2_6)**2) + ((Eq_2_5)**2)); 

Cp :=-.5000000000 (-.2048000000 n 

{.21.0x^-6.x%3j^+15.00x-3.000 + 6.3?%3j^+9.xr^ + 9.3?-3.0i^ + %3j^-3.i^%3jm) 
+ .2048000000 n 

{-6. /^ + 9. 4 ? + 6. jc^ %1 + 9.xr^-3.r^%l + - 6. 4 ? - 6. x %1 ^4? + /^ + %1 + f 


- .06250000000 [-.2048000000 


«(4.00;^ + 4.0;^ + 7.0xH-5.04?-2.x^7^ + 2.x<-6.x% 3?^ j^-4.4-1.0x+3./^%3j^) 

_ n{zx* + i^-2.x^i^+:^ + 3.r^%lJ>^ + r^ -6.x%li^Jj? + r^-4.r^~3.xr^~3.:!^)? t 

+ .2048000000- ^ - f —- I / 

+ ) 

jt^ . .06250000000 (-.2048000000 n 

(-21.0 4 ? - 6. X %3 J%2 + 15.00 X - 3.000 + 6. x2 %3 + 9. X /^ + 9. x3 - 3.0 + %3 - 3.%3 J^) 

lJ^+.204SmOOOn 

(-6.r2 + 9.x3 + 6.4?%l jx? + ;^ +9.x;^-3.r2%l j4? + /^-6.4?-6.x%lJj^ + /2 +%i 

+ 7 ^ y 

%1 := 111 ( 2 . + - 2. x) 

%2:=4?-2.0x+1.00 + 4 
%3:=ln(2.J^ + 2.0-2.x) 
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> with(linalg): 

Warning! new definition for 
Warning: new definition for 


norm 

trace 


> 

>points:= 10; 

> # points:: 48; 

> 

> increment:: (1.0) / (points +1); 


points := 10 


increment := .09090909091 


> 

> X(1):: increment; 


X(l):= .09090909091 


> r(1);= 2.0 * tau * ((X(1)) - (X(1)*2)); 

r(l):= .02644628100 

> nonjdim_ux(1):= evalf( subs( x:X(1), r:r(1), Eq_2_5)); 

non_dim_ux(l) := -.08235051778 

> nonjlim_ur(1):= evalf( 6ubs( x:X(1), r:r(1), Eq_2_6)); 

non_dim_ur( 1) := .2511767847 

> cp(1):: evalf( 8ubs( x:X(1), r:r(1), Cp )); 

cp(l):= .09482965073 

> for i from 2 by 1 to points do 

> 

> X0):= X(i-1) increment: 

> r0):= 2.0*tau*((X(i))-(X(i)*2)): 

> non_dim_ux(i):: evalf( 8Ubs( x:X(i), r:r(i), Eq_2_S)); 

> non~dim_ur®:: evalf( subs( x:X0), r:rO), Eq_2_6)); 

> cpO):: ~ evalf( subs( x:X0), r:rO), Cp 7); 

> 

>od: 


>#writeto(cp); 

> 

> for i from 0 by 1 to points do 

> 

> if (i = 0) then 

> printf(‘ ‘):iprint0: 

> printf( ‘ tau : %4.2f, tau ):iprintO: 

> printf( ‘ ‘ ):iprint0: 

> printf(‘ X r ux ur Cp ‘ ):iprintO: 

> printfC ‘):iprint0: 

> eise 

> printf(‘ %8.4f %8.4f %8.4f %8.4f %8.4f X(0,r(i),non_dim ux(i),non dim urCi),cp(i)):lprintO 

> fi 

> 

>od; 
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tau 


0*16 


X 

r 

ux 

ur 

Cp 

0.0909 

0.0264 

-0.0824 

0.2512 

0.0948 

0.1818 

0.0476 

-0.0136 

0.1874 

-0.0082 

0.2727 

0.0635 

0.0244 

0.1300 

-0.0663 

0.3636 

0.0740 

0.0458 

0.0766 

-0.0996 

0.4545 

0.0793 

0.0557 

0.0253 

-0.1150 

0.5455 

0.0793 

0.0557 

-0.0253 

-0.1150 

0.6364 

0.0740 

0.0458 

-0.0766 

-0.0996 

0.7273 

0.0635 

0.0244 

-0.1300 

-0.0663 

0.8182 

0.0476 

-0.0136 

-0.1874 

-0.0082 

0.9091 

0.0264 

-0.0824 

-0.2512 

0.0948 


># SECTION 4. Calculate the approximate non<dimensional, longitudinal velocity component, 

> # using Equations 2.10 and compare to Equation 2.5. 

> 

>Fp1x:= eq_2_3; 

Fplx := 2.0 (.3201 - .320(.320 - .640 5) 

>Fp2x:= diff( eq_2_3, xl); 

Fp2x := 2.0 Tc (.320 - .640 -1.2800 n (^201 - .320 f) 

> Fp3x:= dlff( eq_2_3, xi$2); 


Fp3x := -3.8400 n (.320 - .6401) 

>Fp4x:= diff( eq_2_3, xi$3); 

Fji>4j::= 2.4576000 Ji 

> tempi := sqrt( xi * (1.0-xl)); 

tempi 

> temp2:= (Fplx) * ((O-S-xl) / (xl*(1.0-xi))); 

temp2 -2.0 "U20|-320|^)(.320..640|)(J.^) 

I(i-o-l) 

> temp3:= (Fp2x) * (in( 2.onemp1) -1.0); 

temp3 ;= ( 2 .Oji (.320 - .640 -1.2800 n (.3201 - .3201^)) (ln( 2.0 J%JT^). 1 . 0 ) 

> temp4:= (Fp3x/ 2.0) * ( 0.5-xl); 


temp4 := -1.920000000 n (.320 - .640 5) (.5 -1) 

> temp5:= (Fp4x/48.0) * (1.0 + (4.0*((0.5-xl)**2))); 

tenq>5 := .05119999999 Jt (1.0 + 4.0 (.5 - 

* E*L,2_1 0:= (Fp2x/(2.0*pl))*(ln(r))-((1.0/(2.0*pi))*(temp2-ftemp3+temp4-i-temp5)); 
EqJJO .5000000000 1.2800. (j20S-320|^)) Ipfr) 

Jl 

^ ^ n (.320 g - .320 1) (-320 - .6401) (.5 - 1) 

1 (10-1) 


-.5000000000 


( 

[ 
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+ (2.0 n (.320 - .640 if ■ 1.2800 n (.320 ^ - .320(ln(2.0 Ji" jl.O 4) - 1-0) 

-1.920000000n (.320 - .640(.5 -1) + .05119999999 jt (1.0 + 4.0 (.5 - )]/ji 


> 

>points;= 10; 

points := 10 

> if points::: 48; 

> 

> increment:= (1.0) / (points + 1); 

increment := .09090909091 

> 

> X(1):= Increment; 

X(l) ;= .09090909091 

> r{1):= 2.0 * tau * ((X(1)) - (X(1)*2)); 

r(l) := .02644628100 

> non_dim_ux(1):= evaif( subs( 3fc:X(1), r=r(1), Eq_2_5)); 

non_dim_ux(l) := -.08235051778 

> approx_eq_2_5(1):= evalf( subs( xi=X(1), r=r(1), Eq_2_10)); 

approx_eq_2_5(l) := -.0815224471 

> 

> for i from 2 by 1 to points do 

> 

> X(i):= X(i-1) + increment: 

> r0):= 2.0*taii*((X(l))-(XG)*2)): 

> non_dim_ux(i):= evalf( 8ubs( x=X0), r=r®, Eq_2_5)); 

> apprbx_eq_2_SG):= evalf( subs( xiaX(i), r=r0), Eq_2_10 )); 

> 

>od: 

> 

> # writeto( eq36ux); 

> 

> for i from 0 by 1 to points do 

> 

. > if G = 0) then 

> printf(‘ ‘):iprint0: 

> printf( ‘ tau = %4.2r, tau ):iprintO: 

> printf(‘ ‘}:lprint0: 

> printf(‘ X r Eq_2_10 Eq_2 5 ‘):ipi1nt0: 

> printfC ‘):iprint0: 

> eise 

> printfC %8.4f %8.4f %8.4f %8.4f,XG),rG),approx_eq_2_5G),non_dim_uxG)):lprint0; 

> fi 

> 

>od; 

tau = 0.16 

X r Eq_2_10 Eq_2_5 
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0.0909 

0.0264 

0.1818 

0.0476 

0.2727 

0.0635 

0.3636 

0.0740 

0.4545 

0.0793 

0.5455 

0.0793 

0.6364 

0.0740 

0.7273 

0.0635 

0.8182 

0.0476 

0.9091 

0.0264 


-0.0815 

-0.0824 

-0.0141 

-0.0136 

0.0222 

0.0244 

0.0423 

0.0458 

0.0514 

0.0557 

0.0514 

0.0557 

0.0423 

0.0458 

0.0222 

0.0244 

-0.0141 

-0.0136 

-0.0815 

-0.0824 


> 


> 

> # SECTION 5. Verify that as r goes to zero and x goes to zero, 

> # Cp values will become Indeterminate. 


> 

> points:= 60; 


points 50 


> 

> increments (1.0) / (10.0 + 1.0 ); 

increment ;= .09090909091 


>X(1):= 


>r(1):= 

>cp(1):= 


Increment; 


X(l) 

2.0*tau*((X(1)).(X(1)"2)); 


r(l): 

evalf( subs( x=X(1), r-r(1), Cp)); 


cp(l) 


.09090909091 

.02644628100 

.09482965073 


> for i from 2 by 1 to points do 

> 

> X(|):= XO-1) * Increment: 

> rfi):= 2.0*tau* ((X(l)) - (X(l)*2)): 

> cpO):= evalf( subs( xsXO), r=r(i), Cp)); 

> 

>od: 

> # writeto( xO); 

> for i from 0 by 1 to points do 

> 

> if (i = 0) then 

> printf( ‘ ‘ ):lprintO: 

> printf( ‘ tau = %4.2f, tau ):lprintO: 

> printf( ‘ ‘ ):lprint0: 

> printf(‘ X r Cp ‘ ):lprint0: 

> printf( ‘ ‘ ):lprintO: 

> else 

> printf(‘%18.16e %18.16e %19.16e‘,X(i),r(i),cp(i)):lprintO; 

> fi 


> 

> od; 


tau = 0.16 
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X 


r 


Cp 


9 .0909090909999990e-02 

8.2644628099999978e-03 

7.51314800899999856-04 

6.83013455400000106-05 

6.20921323099999926-06 

5.64473930100000016-07 

5.13158118300000056-08 

4.66507380299999916-09 

4.24097618500000046-10 

3.85543289499999986-11 

3.50493899500000026-12 

3.18630817700000036-13 

2.89664379700000046-14 

2.63331254299999996-15 

2.39392049399999996-16 

2.17629135800000016-17 

1.97844668900000026-18 

1.79858789899999986-19 

1.63507990800000016-20 

1.48643628000000006-21 

1.35130570900000026-22 

1.22845973500000016-23 

1.11678157700000016-24 

1.01525597900000026-25 

9.22959980899999996-27 

8.39054528099999966-28 

7.62776843700000066-29 

6.93433494300000016-30 

6.30394085699999996-31 

5.73085532500000026-32 

5.20986847700000036-33 

4.73624407000000036-34 

4.30567642700000006-35 

3.91425129700000006-36 

3.55841026999999946-37 

3.23491842699999936-38 

2.94083493400000036-39 

2.67348630400000016-40 

2.43044209499999966-41 

2.20949281400000026-42 

2.00862983100000006-43 

1.82602711899999996-44 

1.66002465400000006-45 

1.50911332200000006-46 

1.37192120199999976-47 

1.24720109300000006-48 

1.13381917499999986-49 

1.03074470499999996-50 

9.37040640900000066-52 

8.51855128100000146-53 


2.64462809999999996-02 

2.62277166800000026-03 

2.40240104600000016-04 

2.18549377599999986-05 

1.98693589700000016-06 

1.80631555699999996-07 

1.64210589399999996-08 

1.49282361000000006-09 

1.35711237899999986-10 

1.23373852600000006-11 

1.12158047800000016-12 

1.01961861699999986-13 

9.26926015000000186-15 

8.42660013800000046-16 

7.66054558000000106-17 

6.96413234600000096-18 

6.33102940399999966-19 

5.75548127599999956-20 

5.23225570599999926-21 

4.75659609600000076-22 

4.32417826800000076-23 

3.93107115199999926-24 

3.57370104600000056-25 

3.24881913199999906-26 

2.95347193800000036-27 

2.68497449000000056-28 

2.44088590000000016-29 

2.21898718200000046-30 

2.01726107400000026-31 

1.83387370400000036-32 

1.66715791300000006-33 

1.51559810200000026-34 

1.37781645700000016-35 

1.25256041500000026-36 

1.13869128599999986-37 

1.03517389700000026-38 

9.41067178800000076-40 

8.55515617199999916-41 

7.77741470400000006-42 

7.07037700399999996-43 

6.42761546000000076-44 

5.84328677999999956-45 

5.31207889199999986-46 

4.82916262999999986-47 

4.39014784599999996-48 

3.99104349800000076-49 

3.62822135999999936-50 

3.29838305600000016-51 

2.99853005000000036-52 

2.72593640999999946-53 


9.48296507300000206-02 

3.61368535200000016-01 

5.35217608200000086-01 

6.67988135099999996-01 

7.69335666599999976-01 

8.40414870299999956-01 

8.81338185900000026-01 

8.92114847999999966-01 

8.72745451000000016-01 

8.23230012999999986-01 

7.43568530999999986-01 

6.33761003000000166-01 

4.93807431999999946-01 

3.23707815999999986-01 

1.23462154999999996-01 

-1.06929549999999976-01 

-3.67467300999999976-01 

-6.58151099000000026-01 

-9.78980939000000166-01 

-1.32995682300000016+00 

-1.71107875800000026+00 

-2.12234672700000056+00 

-2.56376076000000016+00 

-3.03532080799999986+00 

-3.53702692800000006+00 

-4.06887908700000006+00 

-4.63087728199999936+00 

-5.22302152700000026+00 

-5.84531181199999946+00 

-6.49774815200000076+00 

-7.18033052700000066+00 

-7.89305895200000056+00 

-8.63593341700000046+00 

-9.40895394700000056+00 

-1.02121204800000016+01 

-1.10454331000000006+01 

-1.19088917399999986+01 

-1.28024964300000016+01 

-1.37262471900000006+01 

-1.46801439600000016+01 

-1.56641867600000026+01 

-1.66783756199999976+01 

-1.77227105899999986+01 

-1.87971915399999966+01 

-1.99018185400000046+01 

-2.10365915700000026+01 

-2.22015106799999956+01 

-2.33965758100000056+01 

-2.46217870199999996+01 

-2.58771442200000026+01 


> 


> 

> # SECTION 6. Verify that as r goes to zero and x goes to one. 
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Cp values will become indeterminate. 


>X(1):= 0.9500; 

X(l):=.9500 

> r(1):= 2.0 * tau * ((X{1)) - (X(1)'‘2)); 

r(l):= .01520000000 

> cp(1):= evalf( subs( x=X(1), r=r(1), Cp)); 

cp(l) := .1752243365 

>X(2):= 0.999900; 

X(2):= .999900 

> r(2):= , 2.0*tau * ((X(1)) - (X(1)'‘2)); 

r(2) := .01520000000 

> cp(2):= €valf( subs( x=X(1), r=!r(1), Cp)); 

cp(2) := .1752243365 


> points:= 29; 


> incrementrs 0.000001; 


> for i from 3 by 1 to points do 


points := 29 


increment := .110"® 


> X0):s X(i-1) + increment: 

> r(i):= 2.0nau* ((X(i)) - (X(i)'‘2)): 

> cp(i}:= evalf(subs(x=XO), r=r(i),Cp)); 

> 

> od: 

Error, (in evalf/ln) singularity encountered 

>#writeto(x1); 

> for i from 0 by 1 to points do 

> 

> if (i 0) then 

> printf( ‘ ‘ ):lprint0: 

> printf( ‘ tau = %4.2f, tau ):lprintO: 

> printf( ‘ * ):iprint0: 

> printf(‘ X r Cp ‘ ):lprintO: 

> printf( ‘ ‘ ):lprintO: 

> else 

> printfC %18.16f %18.16f %19.16f‘,X(i),r0),cp(0):lprint0; 

> fi 


tau = 0.16 


0.9499999999999998 0.0152000000000000 0.1752243365000000 
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0.9998999999999999 

0.9999009999999999 

0.9999020000000001 

0.9999029999999999 

0.9999040000000001 

0.9999050000000002 

0.9999060000000001 

0.9999069999999999 

0.9999079999999999 

0.9999090000000000 

0.9999100000000000 

0.9999110000000002 

0.9999120000000001 

0.9999129999999999 

0.9999140000000001 

0.9999150000000002 

0.9999160000000001 

0.9999169999999999 

0.9999179999999999 

0.9999189999999999 

0.9999200000000000 

0.9999210000000001 

0.9999219999999999 

0.9999230000000002 

0.9999240000000001 

0.9999250000000000 

0.9999260000000000 

0.9999269999999999 


> # file s app_a.txt 

> 

> 

> 


0.0152000000000000 

0.0000316768640000 

0.0000313569280000 

0.0000310369920000 

0.0000307170560000 

0.0000303971200000 

0.0000300771840000 

0.0000297572480000 

0.0000294372800000 

0.0000291173440000 

0.0000287974080000 

0.0000284774720000 

0.0000281575360000 

0.0000278375680000 

0.0000275176320000 

0.0000271976960000 

0.0000268777280000 

0.0000265577920000 

0.0000262378560000 

0.0000259178880000 

0.0000255979520000 

0.0000252780160000 

0.0000249580480000 

0.0000246381120000 

0.0000243181440000 

0.0000239982080000 

0.0000236782400000 

0.0000233583040000 


0.1752243365000000 

0.6476606865000000 

0.6553688371999999 

0.6422468258999999 

0.6393732819000001 

0.6585344844000000 

0.6561241774000001 

0.6536068717000000 

0.6398706602999999 

0.6489621221999999 

0.6580446579000001 

0.6431658841000000 

0.6400647958000000 

0.6507920520000000 

0.6343043787000000 

0.6448765211000000 

0.6425597706000000 

0.6392634326000001 

0.6357834041000000 

0.6482360795000001 

0.6293168993000000 

0.6416774725000000 

0.6547907977000000 

0.6354180657000001 

0.6494795412000000 

0.6623045039000000 

0.6253918384000000 

5.8453118119999994 
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APPENDIX B. FORTRAN LISTING OF THE SPHERE CODE 



' 0 . 



i it 


IlisI 

iim 

S JZ & 

0<a • f 


c c 
Mje 
c c 

9 9 


II 

*»♦> 

• • 


ii 

m • 


a 


u « ^ 

■ 9 C 

Q.U<H 


if : 

55 e 


S5 

55 

• • 


9 «i 
&U 

• o 


O t6 

• g 

g.g 

ss 

5? 

« b 

10 a 

• a 

55 


4J — 

t 


o 

• g 


d Dt -HU 
w c %4 2 


e 0 0 

« • 

*t c c 

S II 

s 

8 gg 

_ e c 


*i X 
9 Ai 

au 


5? 

« b 

*0 a 




I • 
9 « 


• w 

i 


b ooooooooirtmu^inintn 

9 99999999999999 

9 ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ 

- Mbbbbbbbbbbbbb 


oooooooooooooooooooooo 


U 0 O U O O O 
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0 Caldulftt* thi ■urfao* point*, wh*r*, 17 continue 

o (*) 1* the aphere'* radlua, do 19, j • 1, n 

o (li) la the horiiontal axle of the aphere, dd(n,j) * a(5} 

o jrli la the perpendiouler dlatance freei (tl}« 19 continue 











25^*87 

•.♦TJ , C 

4t wH »- 
9 O 


s • i • u! 

“■x: c 9 

_** c « « 
e - o 


^ .SS£ 

•H O 0*4 4> 
kl «*> O «< ■ 
C C 

0 o o o 

*> T> 0 0 


'"fiS. . 
5*8 


C .-.—cec -c-"- 

O** «*««** ^ -m* g 

^-ooooooe 

0 MW IV <>»« «M 4J 

^ g • « ■ Ji f ■ 0 

.00000000 ^ 

ie 4 >«Jf»vp^v 

^ 0-H-H'►I •H ^ w4<H <M v-k Jif a z 

gUbVibiikiliU « S 

Ow»>>>»»»0 

C 0V ^ Wj< C o 

0V« •-•COOOUOU'-'iaz 

TJ-« 0 0 TJ t3 <0 O V? U K S la 
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function d0t0nn( array, 
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829730996839018R- 
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aph:o«t 


All 20 <99508:39 


S.2249e5647159<89E-O01 

5.877a52522924731E-00l 

€.4944e048330ie37E-001 

7.0710678X1885476E-001 

7.6040S96S6000309E-001 

8.09016994374947SE-001 

8.S26401643540922E-001 

8.910065241883679E-001 

9.238795325112867E-.001 

9.S1056516295X535E-001 

9.723699203976766E-001 

9.876e8340595l378E-001 

9.969173337331280E-001 


8.526401643540922E>001 
8.090I69943749475E->00X 
7.604059658000309E-OOX 
7.07ie67eX18€547SE-00X 
6.494480483301837E»001 
5.e7785252292473XE-001 
S.2249e56471S9488E~001 
4.S399O4997395468E~00I 
3.826834323650898B>-00X 
3.090X«9943749474B>OOX 
2.3344S3€38S59054B>00X 
X.564344<50402309E-O01 
7.845909572784494E>002 


CP# 

8.905393839067629E-001 
9.027353434792928E-001 
e.5XXX79947«85362E-00X 
7.664369495328547E-00X 
6.564562X48X95853E-001 
5.255XOe32e774930E-OOX 
3.7742605788859S4B-001 
2.X«X3843X10680e6B-OOX 
4.5B08182S9550300B-002 
.1.2922471SX880S2«E>001 
.3.0453294779a596OE-00X 
-4.7S706e22785X425E-001 
.6.364S9940529364«E>001 
-7.8873X3331453S88B-001 
-9.227818362308e64B-00X 
-X.037282781301297 
-1.X2939S086758662 
.X.X98836872040474 
-.1.2379378918X508X 
-1.2SX679365O50C49 
-X.2377X8124165071 
.1.196394403699731 
.1.128723005693080 
.1.036368004554203 
-9.2X60X5679375731E.001 
-7.e7247e831464250B-001 
.6.366X3S399787071B.001 
.4.734060719099040B-001 
.3.OX642648313S478B-O0X 
-1.255SX1S79397279E-001 
5.053388664749303B.002 
2.222780643714663E.001 
3.854537S60597832E-001 
S.3604422445SX401E-001 
6.703425006687579E.OOI 
7.8S0426542218033E-001 
8.773211116376738E.001 
9.449068531375339B-001 
9.86X3ie536673343E-001 


cpex 

9.861493831695299E.00I 
9.449385e0832O47?E-OOI 
8.773823397119137B.001 
7.85X441186718157B.001 
6.7049512883486S7E-4>01 
5.362584088290321K-001 
3.8573931220698996.001 
2.226441186718156E-001 
5.098877317025950E.002 
-1.250000000000003B.001 
.3.009807731703601E.OOI 
.4.726441186718159E-001 
.6.3S739312206990SB.001 
.7.8625640S829O32SB.0OX 
-9.2049512883486566.001 
.1.035144118671816 
-1.127382339711914 
.1.194938580832048 
.1.236149383169530 
.1.250000000000000 
.1.236149383169530 
.1.194938580832048 
-1.127382339711913 
-1.035144118671816 
.9.2049512883486566-001 
.7.8625840882903206-001 
.6.3573931220699006-001 
.4.7264411B67I81S9B-4)0X 
.3.0098677317025966-001 
.1.2499999999999986-001 
5.0988773170259726-002 
2.2264411867181586-001 
3.6573931220699d36.001 
5.3625840082903236-001 
6.7049512683486606-001 
7.8514411867181596-001 
8.7738233971191396-001 
9.4493858083204776-001 
9.8614938316952996-001 


In subroutine CRMfER, 
denan» O.OOOOOOOOOOOOOOOB-»-000 
causes {x(k)} in subroutine to blow op, 
CAM MOT obtain a {q) vector at, 

SB 41, su^ace points til 
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APPENDIX D. ERR.OUT FILE FROM THE SPHERE CODE 
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IBND NBC 

$BCIMP lint 


















$SCEIKS> SEND 


$^IDNAM 

HAMB ■ 'Lover Canard Bod • 5, 51x94x33 )•, $BHD 
$KITeRS $BND 
SMETPRK 

IRHS - 0, ILHS - 2, , : 5S - 2, 5EK0 

5TIMACU 

DT «• 0,04, ITIKE- I, CTXKIM-O, CFLKAX«200, TTOSO • X.OO, 523© 
$SM0ACU 

ISPBC - 2, D162 - 2.00, D164 • 0.02, 

SHOO • 1.000, EPSB » 0.0«0, $81© 

$VX8IXP 

VISCJ • .T., VISCK - VISCL » .T., 

NTURB - 2, 

ITWP - 1, 1, 

ITDIR « 3, 1, 

JILS - 1, 1, 

JTLB « SI, 51, 

KTLS « 1, 1, 

KTLE - 94, 94, 

LTLS • I, 1, 

LTLB • 33, 33, 

TLPARl* 0.9, 0.9, 

$END 

SBCXHP 

MBC « 4, 

IBTTP - 15, 10, 5, 5, 

XBDIR • .1, 2, 3, 1, 

JBCS - 51, 1, I, 1, 

JBCB > 51, 51, 51, 1, 

KBCS - 1, 1, 1, 1, 

KBCB - 94, 1, 94, 94, 

1BC8 - 1, 1, 1, 1, 

LBCK - 33, 33, 1, 33, 

581© 

$SCBX1IF $BND 


$ORDIIAK 

WM C ■ 'ORMir Canard Body ( gr_6, 51x94x33 )\ $81© 

$MX1!EH8 $81© 

$»STPRM 

XBB6 • 0; lUS - 2, ZDZS8 « 2, $BHD 
$TZIttCU 

DT ■ 0.04, ITIKB- 1, CnJKZll«0, CPU(AX*200, TF060 ■ 1.00, $81© 
$SM0JkCU 

Z8P8C - 2, DZ82 - 2.00, 0X84 - 0.02, 

8MCX> > 1.000, 8P8B « 0.060, $880 
$VX8X1IP 

VI8CJ - .1., VISCX • .T., VI6C1 .T,, 

WTURB » 2, 

ITTZP * 1, 1, 

ZTDZR « 3, 1, 

2TL8 • 1, 1, 

JT18 - SI, 51, 

KTLS - 1, 1, 

KTL8 « 94, 94, 

LTLS - 1, 1, 

LTLB • 33, 33, 

TLPARl- 0.9, 0.9, 

$81© 

$BCZMP 

HBC - 4, 

ZBTItP - 15, 10, 5, 5, 

XaDZR - -1, 2, 3, 1, 

JBCS « 51, 1, 1, 1, 

JBCB - 51, 51, 51, 1, 

KBCS - 1, 1, 1, 1, 

XBCB - 94, I, 94, 94, 

■ 1, 1, 1, 1, 

LBC8 • 33, 33, 1, 33, 

$BHD 

$SCBZlfP $BND 
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